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ABSTRACT 

Ab-initio numerical study of collisionless shocks in electron-ion unmagnetized plasmas is performed 
with fully relativistic particle in cell simulations. The main properties of the shock are shown, focusing 
on the implications for particle acceleration. Results from previous works with a distinct numerical 
framework are recovered, including the shock structure and the overall acceleration features. Particle 
tracking is then used to analyze in detail the particle dynamics and the acceleration process. We 
observe an energy growth in time that can be reproduced by a Fermi-like mechanism with a reduced 
number of scatterings, in which the time between collisions increases as the particle gains energy, and 
the average acceleration efficiency is not ideal. The in depth analysis of the underlying physics is 
relevant to understand the generation of high energy cosmic rays, the impact on the astrophysical 
shock dynamics, and the consequent emission of radiation. 

Subject headings: acceleration of particles - collisionless shocks - gamma rays: bursts 



1. INTRODUCTION 

The dynamics of relativistic particle acceleration in 
collisionless shocks is of great importance to several 
astrophysical scenarios. The acceleration of electrons, 
positrons, or ions in various structures such as active 
galactic nuclei, gamma-ray bursts, pulsar wind nebu- 
lae, and supernova remnants results in energetic particles 
that can then scatter in magnetic and electric fields and 
emit synchrotron radiation (see, for instance, Jones & El- 
lison 1991). Despite their relevance to understanding the 
radiation collected in astronomical observations, the un- 
derlying processes inherent to the acceleration are not yet 
fully understood. In this context, full kinetic simulations 
play an important role in the assessment of particular 
physical mechanisms relevant for astrophysical shocks, 
namely in t he study of the nonlinea r growth of the W eibel 
instability (|Weibel||1959| fMedvede v fc Loeb||199 9D with 
magnetic field generation (Silva et al. 2003; Fonseca et al. 
2003[[FYederiksen et al. 2004; Nishikawa et al 2 005), and 
particle acceleration in the percursor region (|5ilva 2006 ) . 
Nevertheless, the self-consistent modeling of a relativistic 
collisionless shock, from first principles, is computation- 
ally very demanding, as large temporal and spatial scales 
push the numerical techniques to conditions still unex- 
plored. Hence, full kinetic simulations require massive 
computational resources, optimized algorithms, methods 
for improved energy conservation, and also advanced vi- 
sualization diagnostics. 

Recently, progress was made in the understanding of 
electron-ion shock formation and electron-po sitron accel- 
eration in relativistic unmagnetized shocks fSpitkovsky 
2008a|b ). These studies with particle in cell (PIC) sim- 
ulations confirmed the capability of these structures to 
effectively accelerate electrons, which is identified by the 
development of a non-thermal tail in the energy spec- 
trum. Similar approaches with large-scale self-consistent 
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modelling can provide valuable input to improve Monte 
Carlo methods (e.g., Bednarz & Ostrowski 1998, Ellison 
& Double 2004), and to support the development of an- 
alytical models (e.g.. Kirk et al. 2000, Achterberg et al. 
2001, Keshet 2006). 

Here, we examine a relativistic electron-ion unmagne- 
tized shock with ab-initio relativistic PIC simulations, 
and apply a full particle tracking diagnostic to better 
understand the particle dynamics and the acceleration 
mechanism. We then show that non-thermal particle ac- 
celeration occurs through a small number of scatterings 
in the shock front, in a Fermi-like process (A£^ = aE) 
with increasing time between each energy gain. Our sim- 
ulations and initial data analysis follow and confirm pre- 
vious results obtained with a different PIC framework 



( |Spitkovsky|2008a^b) , which reveals the robustness of the 
numerical results with distinct algorithms and implemen- 
tations. In §2, we present the simulation results of the 
electron-ion plasma shock with a reduced ion to electron 
mass ratio of 32. In §3, we leverage on the OSIRIS ( |Fon-| 
seca et al. 2002) particle tracking and data processmg 
tools to analyze the particle dynamics and their acceler- 
ation process, by focusing on the time evolution of the 
main physical quantities of the most energetic particles. 
A discussion of the acceleration mechanism and the con- 
clusions are presented in §4. 

2. SHOCK FORMATION AND EVOLUTION 

Numerical simulations were performed with OSIRIS, 
a fully relativistic, electromagnetic, and massivelly par- 
allel PIC code which has been used in many different 
physical scenarios, such as astrophysics (e.g., Silva et 
al. 2003), laser/plasma accelerators (e.g.. Mangles et al. 
2004, Tsung et al. 2004), nanoplasma dynamics (e.g., 
Peano et al. 2006), and fast-ignition (e.g., Ren et al. 
2004). 

We simulate a two-dimensional system of a cold 
unmagnetized electron-ion plasma with mass ratio 
mi/rrie = 32 (m^ and rrie the ion and electron mass, 
respectively), and evolve it until evidence of a non- 
thermal acceleration tail in the downstream particle spec- 
trum is achieved. To generate the shock, the plasma 
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stream is launched from the right wall with proper ve- 
locity u = 7/3 = —20, and minimal thermal dispersion 
from randomized particle injection. This neutral plasma 
stream is reflected from a rigid boundary at the left 
wall (this is one of the most direct methods to gener- 
ate shocks in simulations; see, for instance, Forslund et 
al. 1970, Jones & Ellison 1991, or Spitkovsky 2008a,b). 
The computational domain is 50 c/cjp in the transverse 
direction and 280 c/cOp in the longitudinal direction, with 
c/ujp — {^'Ke^nl^mi(?^~^l^ the ion skin depth for a num- 
ber density n and relativistic ion mass 7771^; e is the el- 
ementary charge and c the speed of light in vacuum. A 
time step of 0.012/cc;p is used. The system is numeri- 
cally resolved with 10 cells per electron skin depth in 
both directions, thus ensuring that the dynamics of the 
lighter species is accurately m odeled. Taking a dvantage 
of second order particle shapes (Esirkepov 2001 ) and cur- 
rent smoothing compensation, we use 2 particles per cell 
(ppc) for each species, which is equivalent to 6 ppc in the 
shocked gas, for which the energy conserva tion is equiva- 



lent t o 16 ppc with linear particle shapes ( jFonseca et al. 
20081). Lower mass ratios, higher grid resolutions, and 



more particles per cell were tested, showing an overall 
qualitative and quantitative result convergence. 

The main physical processes that generate the shock 
are observed when the reflected particle stream interacts 
with the incoming plasma, which leads to the growth of 
the Weibel instability, particle thermalization, and the 
generation of electric and magnetic turbulence. The tur- 
bulence slows down the flow, which generates the shock 
as a density compression that propagates in the posi- 
tive x\ direction (Fig.jl^). Simulation results agree with 



the hydro dynamical jump co nditions ( |Blandford fc Mc _ 
Keel 1976 Spitkovsky||2008a ); the steady state velocity 
of tne shock, /:^shock — 0.48, and the corresponding den- 
sity compression obtained directly by particle number 
conservation, 712 /ni = 3.1. We emphasize that, given 
the simulation configuration, all quantities discussed are 
measured in the downstream frame. 

Fig. [1] also shows other relevant physical quantities 
after the shock has achieved a steady state. Of rel- 
evance to the acceleration process, transverse electric 
fields (Fig. [i]d) arise in the linear stage of the Weibel 
instability, associated with space charge effects, as the 
two counter-propagating plasma streams pinch/filament 
with d ifferent rates because of their different tempera- 
tures (Tzoufras et al. 2006). The spatial symmetries 
of this field have direct impact on the overall trans- 
verse momentum acquired by the particles (c.f. 
In addition, the energy deposited in the magnetic neid 
reaches cb = B'^ /ATrjnmc^ 15 — 20%, similarly to 
pair shocks ( Spit kovsky||2008b ) , and to the mi/rrie = 16 
case. This observation, coupled with the average value 
of Vdrift = E X across the shock front associated 

with the structure of the self-consistent fields, suggests 
also the origin of the particle trapping mechanism in the 
shock front for both positrons/ions and electrons (Mar- 
tins et al., in prepa ration) , simila r to what is observed in 
Earth bow shocks (Burgess||2007). 



In accordance with results obtained by Spitkovsky 
(2008b), the energy spectrum of the ions (Fig. [l]^-f) 
is significantly different across the longitudinal direc- 
tion. The upstream region is dominated by the quasi- 
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Fig. 1. — Steady state structure of a mi/rrie = 32 collisionless 
shock dit t = 360 /ujp after first counter-streaming interaction, a) 
Transverse average of the ion density; b) Transverse electric field; 
c) Ion energy spectrum along simulation box; d-f) Ion spectra for 
the downstream, the shock front, and the upstream regions, re- 
spectively. The downstream inset includes also the electron spec- 
trum (scaled down by the mass ratio rrie/mi) with relativistic 
Maxwellian fits for both species - blue dashed lines -, and a fit 
to the electron spectrum of a relativistic Maxwellian with a power 
law trimmed by an exponential cut-off for higher energies - red 
dotted line. 

monoenergetic negative flow of particles, and contains a 
lower density returning stream of heated particles that 
already escaped the shock region (or were never trapped). 
Despite the strong thermalization of the shocked gas, ev- 
idence of non-thermal particle acceleration can be ob- 
served in the downstream ion spectrum (Fig. [l]i), where 
a fit to a pure relativistic Maxwellian does not account 
for the high-energy tail. The non-thermal spectrum of 
both electrons and ions can be fitted with a power law 
(7~^, with p = 2.3 — 2.6) and an exponential cut-off 
defined by exp[-(7 - 7cut)/A7cut], with 7cut = 80 and 
^Tcut = 15 for the ion species. Electrons reach energies 
that are higher by the mass ratio m^/me, thus spanning 
more than three orders of magnitude in energy. The 
high-energy tail (for 7 > 40) has ~2% of the total num- 
ber of ions in the downstream slice analyzed, and ac- 
counts for ~10% of the total ion energy in that spatial 
region. This confirms previous results obtained for a pair 
plasma configuration by Spitkovsky (2008b), leading to 
the conclusion that the generated spectra and acceler- 
ation efficiency are not very sensitive to the mass ratio 
of the species, at least in two-dimensional simulations. 
Finally, additional propagation time leads to a linear in- 
crease of the non-thermal tail span to higher energies, 
and results indicate a fit with the same power law index. 

3. PARTICLE DYNAMICS Sz ACCELERATION 

For a detailed analysis of the particle dynamics and ac- 
celera tion mechanism, the OSIRIS particle tracking fea- 
tures ( Fonseca et al.|20Q8 ) were used to follow the trajec- 
tories of the most energetic particles, selected in a first 
scanning simulation. 

Acceleration occurs for both longitudinal and trans- 
verse momentum, but with different dynamics in each 
direction (Fig. [2|. On one hand, the approximately null 
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space average of the transverse electric field (Fig. [Tp) 
leads to a symmetric acceleration of the particles in the 
X2 direction (Fig. [2]3-d). On the other hand, the longi- 
tudinal momentum shows an average increase over time, 
and the particles reach farther into the upstream region, 
until escaping (Fig. [2k) mostly to the downstream region. 
The propagation oi the particles encompasses several 
magnetic rotations in which momentum is transferred 
between directions, leading to the oscillations observed 
in ui X t and U2 x t. The ions are able to remain inside 
the accelerating structure by performing long drifts in 
the transverse direction {u2 ^ ui). For the most ener- 
getic particle (colored trajectory), the final energy gain 
occurs with high angle variation from the transverse to 
the longitudinal direction (u2 nearly constant in Fig. 
d); the particle is then refiected from the upstream, and 
escapes the shock region into the downstream plasma 
(final ui < 0). After escaping, the particle sustains a 
large transverse momentum, and performs long drifts in 
the transverse direction with constant energy. Also, it 
should be emphasized that the behavior observed for this 
particle is representative of the accelerated particles. 

The energy gain and the interaction with the shock 
region are depicted in Fig. [3^, which shows the time evo- 
lution of the longitudinal position relative to the shock 
front for the 80 most energetic ions. Identically to the 
particle motion in pair plasmas observed in Spitkovsky 
2008b, ions gain energy after being trapped as they per- 
form multiple oscillations in the shock region until they 
finally escape, mainly to the downstream region. The 
wall refiections observed in Fig. [3^ are close to the in- 
jection point only at the beginning of the simulation, 
and do not affect the overall process. After the shock is 
formed, particles from the unshocked gas can be directly 
trapped, without reaching the downstream. The parti- 
cle gains energy from the electric fields of the upstream, 
and then crosses the shock region until being refiected in 
the downstream (Fig.lsb). We emphasize that, since the 
simulation is performea in the downstream frame, no sig- 
nificant acceleration occurs when particles are refiected 
on the downstream shocked gas. Accelerations occur 
rapidly and /S.E :^ £^ is typically observed (Fig. |3|d), as 
expected in a Fermi process (Fermi 1949). The maximum 
energy reached is 7finai — 170. Since for a single bounce 
and the initial ion energy is 7initiai — 18, that fi- 
nal energy can only be achieved after several shock cross- 
ings. It is important to notice that, unlike a Fermi mech- 
anism, the final energy is usually obtained after only 3-5 
effective collisions in a continuously evolving shock. 

In the standard formalism of Fermi acceleration, the 
energy evolution is written as E{N) = £^oexp(A/'), with 
£^0 the initial energy and N the number of energy gains 
(assumed very large: 1), and its time depen- 

dence can therefore be estimated by relating N with 
time t. The initial estimate by Fermi assumed a con- 
stant time between energy gains TcoIi = tq, thereby lead- 
ing to an exponential growth of the energy with time: 
E{t) = £^0 exp(t/ro). To account for the small number of 
discrete energy gains observed in the simulation, we now 
write E{N) = Eo{l^a)^, where we include the constant 
fractional energy gain a < 1 in the form AE = aE^ al- 
ready considered by Fermi. The value a was estimated 
with an individual analysis of the energy gains of each 
tracked particle. Fitting the data of 52 of the 80 tra- 
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Fig. 2. — Momentum-time trajectories of the 80 most energetic 
particles, chosen at t = 360 /ujp and evolved until t = SSO/cjp. 
The trajectory in color corresponds to the most energetic particle, 
a) Longitudinal momentum time evolution (ui x t), b) Longitudi- 
nal/transverse momentum time evolution (ui x U2 x t), c) Longi- 
tudinal/transverse momentum evolution {u2 x ui), d) Transverse 
momentum time evolution (u2 x t). 

jectories to (1 + a)^ yields a o:^ 0.81 (particular energy 
fiuctuations of the remaining particles did not allow for a 
clear identification of all scatterings, defined as xi prop- 
agation direction inversion with > 50% energy gain - see 
arrows in Fig. Isb). 

Further analysis of several particle trajectories in the 
simulation shows an effect that decreases the energy 
growth with time, namely the increase of time between 
energy gains as the particle accelerates and the shock 
structure evolves. A simple model can be obtained us- 
ing Tcoii = To + 5t, for To the initial collision time, and 
where s is the rate of change of the time between colli- 
sions. We thus get an approximate fit N{t) = t/ (tq + 5t), 
and the energy evolution in time becomes E{t) o:^ Eo{l-\- 
Q/)^/('ro+s^) YVe emphasize that this expression is only 
valid for a limited time t, and thus implicitly assumes 
particles escape the accelerating region with a maximum 
of N{t oo) = 1/s collisions. Fig. ^ presents a fit to 
E{t) with s = 0.11 and tq = 50.6/a;p. Alternatively, the 
collision time can be written as a function of energy with 
Tcoii = ^o[E{t)/Eo]^ (an energy dependence of TcoIi is also 
observed in other scenarios, as in the Earth's bow shock 
acceleration, Kis et al. 2004). For this case, a numerical 
fit yields TcoIi oc E{t)^-'^'^^ for the same tq, which assumes 
no time domain restrictions, as opposed to TcoIi oc t. The 
parameters of the model, particularly the fractional en- 
ergy gain a, and the growth rate s, are very similar to 
those obtained with a mass ratio of 16. Nevertheless, the 
parameter study required to explore these dependencies, 
and to completely understand the microphysics under- 
lying the parameter values, is beyond the scope of the 
present paper and will be tackled in future work. 

4. CONCLUSIONS & DISCUSSION 

Ab-initio full PIC simulations have been presented for 
a two-dimensional relativistic collisionless shock propa- 
gating in an initially unmagnetized electron-ion plasma 
(mass ratio of 32). The shock structure and jump con- 
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Fig. 3. — Ion dynamics and energy evolution in time, a) Longi- 
tudinal ion position relative to the shock front (xs = xi — Xghock) 
as a function of time (color in trajectories associated with energy) 
for the 80 most energetic particles (chosen at t = 360/a;p), b) En- 
ergy as a function of the position relative to the shock front for the 
most energetic particle; arrows indicate energy gains, c) Energy 
evolution in time for the same 80 ions (gray) with t = the trap- 
ping moment of each particle, corresponding average (solid line), 
and fit of E{t) ~ Eo{l + a)*/(^o+5t) with a = 0.81, s = 0.11, and 
To = bO.G/ujp (dashed line). The final fiattening occurs as particles 
escape the shock, and leave the accelerating region. 

ditions confirm previ ous results obtained with a differ- 
ent PIC framework (Spitkovsky 2QQ8a|b ). Non-thermal 
particle acceleration is observed and occurs as particles 
are trapped and oscillate in the shock front, similarly 
to a Fermi acceleration process. Nevertheless, specific 
distinctions exist from the standard Fermi mechanism, 
namely the small number of scatterings and the contin- 
uous evolution of the shock structure where the particle 
accelerates. When gaining energy, the particle usually 
performs a rotation to the transverse direction, and thus 
remains in the shock region, being susceptible to further 
acceleration. 

An important consequence of the overall acceleration 
efficiency (~ 10% of energy carried by the most ener- 
getic particles) is the increased relevance of nonlinear 
effects for the shock structure and for the acceleration 
process. In fact, when the accelerated particles yield 
> 10% of the plasma energy, their dynamics becomes rel- 
evant and influences the evolution of the overall system 



(Jones & Ellison 1991). One of these nonlinear effects 
is the dynamic pressure of the accelerated particles that 
slows down the unshocked plasma before it reaches the 
sharp shock transition. This is indicated in the simula- 
tion results as the shock widens and the magnetic and 
electr ic field layers extend to the upst ream (Keshet et al.| 
2008[ [Medvedev fc Zakutnyaya||20Q8| . Furthermore, the 
inclusion of the trapped particles population can have 
implications on the jump conditions, and current mod- 
els can be extended to incorporate these effects into the 
equation of state in the shock region. This generalization 
has been m ade for the non-relativist ic and elec trostatic 
shocks case dForslund et"ar]|1970[ [Sorasio et al.| f2006 ) . 

The average energy growth of tne particles can be re- 
produced by a multi-scattering acceleration mechanism 
with AE = aE^ assuming an increase of the time be- 
tween collisions. The long term acceleration implica- 
tions of these effects cannot be inferred from our sim- 
ulation spectra because of the exponential cut-off due to 
the finite simulated size that limits the injection of par- 
ticles. Larger scale simulations with longer propagation 
distances and larger acceleration times will elucidate the 
properties of the particle spectrum at higher energies, 
thus allowing for a detailed identification of the mecha- 
nisms responsible for a < 1, and for an increase of time 
between collisions. 

In summary, our results confirm the possibility of par- 
ticle acceleration through a Fermi-like mechanism with a 
reduced number of energy gains, and generalized to re- 
produce the statistical data obtained with ab-initio full 
PIC simulations that self-consistently resolve the turbu- 
lent and non-linear evolution of the shock. 
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